function X = fft_custom(x)
    N = length(x);
    % Calculate the smallest power of 2 greater than or equal to N
    M = 2^(nextpow2(N));
    x_padded = [x, zeros(1, M - N)];  % Zero-pad the input signal
    if M <= 1
        X = x_padded;
        return;
    end
    % Recursively compute FFT for even and odd parts
    X_even = fft_custom(x_padded(1:2:end));
    X_odd = fft_custom(x_padded(2:2:end));
    % Precompute the twiddle factors (complex exponential terms)
    T = exp(-1i * 2 * pi * (0:(M/2-1)) / M);
    % Combine even and odd parts
    X = [X_even + T .* X_odd, X_even - T .* X_odd];
end

%%%%%%%%图形化界面%%%%%%%%
function fft_gui()
    % Create a figure window
    hFig = figure('Name', 'FFT Custom GUI', 'NumberTitle', 'off', 'Position', [100, 100, 600, 400]);
    % Input signal text box
    uicontrol('Style', 'text', 'Position', [20, 340, 100, 30], 'String', 'Input Signal:');
    hInput = uicontrol('Style', 'edit', 'Position', [130, 340, 300, 30]);
    % Button to compute FFT
    uicontrol('Style', 'pushbutton', 'Position', [450, 340, 100, 30], 'String', 'Compute FFT', ...
              'Callback', @(src, event) compute_fft(hInput));
    % Axes for plotting
    hAxes = axes('Position', [0.1, 0.1, 0.8, 0.6]);
    title(hAxes, 'FFT Result');
    xlabel(hAxes, 'Frequency (Hz)');
    ylabel(hAxes, 'Magnitude');
end

function compute_fft(hInput)
    % Get the input signal from the text box
    inputStr = get(hInput, 'String');
    inputSignal = str2num(inputStr); %#ok<ST2NM> % Convert string to numeric array
    % Check if input is valid
    if isempty(inputSignal)
        errordlg('Invalid input! Please enter a numeric array.', 'Input Error');
        return;
    end
    % Compute FFT using the custom function
    X = fft_custom(inputSignal);
    % Compute frequency axis
    N = length(X);
    f = (0:N-1) / N;  % Normalized frequency (0 to 1)
    % Plot the magnitude of the FFT
    axes(findobj('Type', 'axes')); % Focus on the correct axes
    plot(f, abs(X));
    title('Magnitude of FFT');
    xlabel('Normalized Frequency (cycles/sample)');
    ylabel('Magnitude');
    grid on;
end

% % 图形化界面
% fft_gui
%%%%%%%%%%%%%%%%%%%%%%%%%

% %%%%%%%加载mat文件%%%%%%%%
% % 加载一个mat数据 这里使用的文件来自于西储大学轴承数据集
% load('B014_0.mat');
% x = X185_DE_time;
% x = transpose(x); % 将x转换为行向量
% %%%%%%%%%%%%%%%%%%%%%%%%%%
% 
%%%%%%%%%生成信号%%%%%%%%%%
% % 参数设置
% Fs = 1000;       % 采样频率（Hz）
% f1 = 50;         % 第一个正弦波频率（Hz）
% f2 = 150;        % 第二个正弦波频率（Hz）
% f3 = 300;        % 第三个正弦波频率（Hz）
% T = 1;           % 信号持续时间（秒）
% t = 0:1/Fs:T-1/Fs; % 时间向量，注意这里减去1/Fs以确保T秒
% 
% % 生成三个正弦波
% x1 = sin(2 * pi * f1 * t);  % 第一个频率的正弦波
% x2 = sin(2 * pi * f2 * t);  % 第二个频率的正弦波
% x3 = sin(2 * pi * f3 * t);  % 第三个频率的正弦波
% 
% % 将三个正弦波叠加
% x = x1 + x2 + x3;
%%%%%%%%%%%%%%%%%%%%%%%%%%

% 使用自定义FFT函数
X_custom = fft_custom(x);
K = length(X_custom);
N = 2^(nextpow2(K));
% 使用MATLAB内置的FFT函数
X_builtin = fft(x,N);
% 计算频率轴
f_axis = (0:N-1) * (Fs/N);


plot(f_axis, abs(X_custom), 'b');
title('Magnitude of Custom FFT');
xlabel('Normalized Frequency (cycles/sample)');
ylabel('Magnitude');
grid on;
% Plot the magnitude of the Built-in FFT
figure;
plot(f_axis, abs(X_builtin), 'r');
title('Magnitude of Built-in FFT');
xlabel('Normalized Frequency (cycles/sample)');
ylabel('Magnitude');
grid on;
